#ifndef geometry_H
#define geometry_H

#include <vector>
#include <bitset>
#include <math.h>

#include "excepts.h"

namespace ABAN
{
    typedef double FloatType;
    typedef int SeqType;

    const FloatType zeroRadious = 1e-6;

    template <class T> inline T sqr(const T param){return param*param;}

    const SeqType noOfAxis = 3;
    const SeqType noOfDirections = 6;

    //Cartesian Directions
    const SeqType dirXP = 0;
    const SeqType dirYP = 1;
    const SeqType dirZP = 2;
    const SeqType dirXM = dirXP + noOfAxis;
    const SeqType dirYM = dirYP + noOfAxis;
    const SeqType dirZM = dirZP + noOfAxis;

    const SeqType oppositeSide[noOfDirections] =
      {dirXM,dirYM,dirZM,dirXP,dirYP,dirZP};

    //Cartesian Mesh Constants
    const SeqType nodesPerCell = 8;
    const SeqType facesPerCell = noOfDirections;
    const SeqType nodesPerFace = 4;

    //Base Vector Structure
    struct Vector;
    typedef std::vector<Vector> VectorContainer;
    typedef VectorContainer::iterator VectorIterator;

    struct Vector
    {
        private:
        FloatType x,y,z;

        public:
        Vector():x(0),y(0),z(0){}
        Vector(FloatType ax, FloatType ay, FloatType az):x(ax),y(ay),z(az){}

        inline FloatType&operator[](SeqType dir);
        bool operator==(Vector &otherVec);
        void operator+=(Vector &otherVec);
        void operator*=(FloatType r);

        static bool compare(VectorIterator a,VectorIterator b);
        static FloatType dist(Vector &aVec, Vector &bVec);
    };

    const SeqType nodeNumbersInFaces[facesPerCell][nodesPerFace]=
    {
        {0,2,3,1},//dirXP
        {0,4,5,1},//dirYP
        {0,2,4,6},//dirZP
        {4,6,7,5},//dirXM
        {2,6,7,3},//dirYM
        {1,3,5,7} //dirZM
    };

    const SeqType facesOfNodes[nodesPerCell][noOfAxis]=
    {
        {dirXP,dirYP,dirZP},  //0
        {dirXP,dirYP,dirZM},  //1
        {dirXP,dirYM,dirZP},  //2
        {dirXP,dirYM,dirZM},  //3
        {dirXM,dirYP,dirZP},  //4
        {dirXM,dirYP,dirZM},  //5
        {dirXM,dirYM,dirZP},  //6
        {dirXM,dirYM,dirZM}   //7
    };
    //a Container of Vector iterators, containing nodes of each cell
    //this structure will be used only in finding matching cell neighbors
    typedef std::vector<VectorIterator> RawVectorIteratorContainer;
    class CellConnectivity : public RawVectorIteratorContainer
    {
        public:
        std::bitset<facesPerCell> flag;

        CellConnectivity():RawVectorIteratorContainer(nodesPerCell),flag(0){}; //default constructor
        void rearrangeNodes();
        bool isMatchingNeighborOnFace(CellConnectivity aGCell, SeqType face);
    };

    //Fractional coordination component
    enum FCRelation{fcrNone,fcrMajor,fcrMinor,fcrEqual};
    struct FractionalCoordinationComponent
    {
        SeqType numerator;//surat :p
        SeqType denominator;//makhraj :p

        FractionalCoordinationComponent():numerator(1),denominator(1){}
        FractionalCoordinationComponent(SeqType num, SeqType denom):numerator(num),denominator(denom){}
        void split(bool firstSec=true);
        void collapse();
        FCRelation relation(FractionalCoordinationComponent aFC);
    };
    typedef FractionalCoordinationComponent FractionalCoordination[noOfAxis];

    /*******************************************************/
    /*************************Vector************************/
    /*******************************************************/
    FloatType& Vector::operator[](SeqType dir)
    {
        return *(&(this->x)+dir);
    }
    //-------------------------------------------------------
    bool Vector::compare(VectorIterator a,VectorIterator b)
    {
        if (fabs((*a)[dirXP] - (*b)[dirXP]) > zeroRadious) return (*a)[dirXP] > (*b)[dirXP]; else
        if (fabs((*a)[dirYP] - (*b)[dirYP]) > zeroRadious) return (*a)[dirYP] > (*b)[dirYP]; else
        if (fabs((*a)[dirZP] - (*b)[dirZP]) > zeroRadious) return (*a)[dirZP] > (*b)[dirZP]; else
            throw new EInputErr();
    }
    //-------------------------------------------------------
    FloatType Vector::dist(Vector &aVec,Vector &bVec)
    {
        return sqrt(sqr(aVec.x-bVec.x)+sqr(aVec.y-bVec.y)+sqr(aVec.z-bVec.z));
    }
    //-------------------------------------------------------
    bool Vector::operator==(Vector &otherVec)
    {
        return dist(*this, otherVec)<zeroRadious;
    }
    //-------------------------------------------------------
    void Vector::operator+=(Vector &otherVec)
    {
        for (int i=dirXP; i<noOfAxis; i++)
            operator[](i)+=otherVec[i];
    }
    //-------------------------------------------------------
    void Vector::operator*=(FloatType r)
    {
        for (int i=dirXP; i<noOfAxis; i++)
            operator[](i)*=r;
    }
    /*******************************************************/
    /*****************FractionalCoordination****************/
    /*******************************************************/
    void FractionalCoordinationComponent::split(bool firstSec)
    {
        denominator*=2;
        numerator*=2;
        if(firstSec) numerator--;
    };
    //-------------------------------------------------------
    void FractionalCoordinationComponent::collapse()
    {
        denominator/=2;
        if (numerator%2) numerator++;
        numerator/=2;
    };
    //-------------------------------------------------------
    FCRelation FractionalCoordinationComponent::relation(FractionalCoordinationComponent aFC)
    {
        if ((numerator == aFC.numerator)&&(denominator == aFC.denominator)) return fcrEqual;
        if /*(numerator != aFC.numerator)&&*/(denominator == aFC.denominator) return fcrNone;

        /* if ((numerator != aFC.numerator)&&(denominator != aFC.denominator)) */
        bool iAmMajor = (denominator<aFC.denominator);
        FractionalCoordinationComponent maj,min;
        if (iAmMajor) {maj=*this;min=aFC;} else {maj=aFC;min=*this;}
        while (maj.denominator != min.denominator) min.collapse();
        return (maj.numerator==min.numerator?(iAmMajor?fcrMajor:fcrMinor):fcrNone);
    }
}
#endif
